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Abstract 

We characterize the radial migration of stars in the disk plane by calculating the diffu- 
sion coefficient and the diffusion time-scale for a bulge-disk A/^-body self-consistent system 
with a marginally- stable Toomre-Q parameter. We find that diffusion is not constant in time, 
but follows the evolution of the bar, and becomes maximum near the corotation region and 
in the external disk region, where asymmetric patterns develop. 

Keywords Galaxies: kinematics and dynamics. Galaxies: stellar content. Galaxies: spiral. 
Galaxy: disk 

1 Introduction 

The majority of local galaxies are barred-spiral galaxies, such as our Milky Way |[T1 O 
where the bar is the strongest non-axisymmetric pattern in the disk. The bar is generally 
coupled to other non-axisymmetric patterns such as spirals and warps if the disk is suffi- 
ciently cold lO. The bar has a time-dependent activity, with a pattern speed which typically 
decreases in isolated galaxies. However, the system can be cooled again by adding dissipa- 
tive infall of gas, or by forming stars on low-velocity dispersion orbits, with the net effect 
of restoring the amplitude of spiral waves and the strength of the bar, or even destroying it. 
In this way bars (and spiral waves) can be seen as recurrent patterns which can be rebuilt 
during their long history until the present configuration at z = [4J. 

Under the action of these non-axisymmetric patterns, stars move in the disk which 
gradually becomes hotter. Velocity dispersion of disk stars rises with age, as confirmed 
by observations in the Solar neighborhood (see [5| and references therein) and in external 
galaxies [i6jl71. The origin and the amount of disk heating are still open to debate. 

First attempts to explain such a heating process in the disk of galaxies tried to model 
empirically the observed increase of the stellar velocity dispersion with age in the solar 
neighborhood. Wielen lO suggested a diffusion mechanism in velocity space, which gives 
rise to typical relaxation times for young disk stars of the order of the period of revolution 
and to a deviation of stellar positions of 1.5 kpc in 200 Myr. The result was obtained with- 
out making detailed assumptions on the underlying local acceleration process responsible 
for the diffusion of stellar orbits. Global acceleration processes, such as the gravitational 
field of stationary density waves or of central bars with constant pattern speed, were ruled 
out since their contribution to the velocity dispersion of old stars is negligible and concen- 
trated in particular resonance regions [8, 9|. Different local accelerating mechanisms have 
been investigated so far in isolated galaxies, such as the gravitational encounters between 
stars and giant molecular clouds |T0l[TTl[l2|, secular heating produced by transient spiral 
arms ||T3l[l4l[T5l or the combination of the two processes l T6HT7ll . 

The radial excursion predicted by Wielen [ 8 1 is not sufficient to explain the weakness 
of the correlation between age and metallicity in the Solar neighborhood (see, for example, 
ifTSI ). In order to explain both the large scatter in the age-metallicity relationship and the 
evidence that even old disk stars today have nearly circular orbits, Sellwood & Binney fW\ 
have recently suggested a new mechanism based on the resonant scattering of stars under 
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the effect of transient spiral waves. In this process, a star initially on a nearly circular 
orbit resonates with a rotating wave and changes its angular momentum. If the duration 
of the peak amplitude of the perturbing potential is less than the period of the 'horseshoe' 
orbits, i.e. orbits of particles trapped at the corotation radius of the spiral wave, the star 
can escape from the potential well without changing its eccentricity. The net effect of this 
scattering mechanism is that stars migrate radially without heating the disk. In other words, 
the overall distribution of angular momentum is preserved, except near the corotation region 
of the transient spiral wave, where stars can have large changes of their angular momenta. 

The mechanisms driving radial migration and heating are still hotly debated. Minchev et 
al. EOlllTI propose a mixing mechanism where resonances between the bar and the spiral 
arms can act much more efficiently than transient spiral structure, dramatically reducing the 
predicted mixing timescales. 

Most of the observational signatures of radial mixing reported in the literature f22l [23l 
[24ll point to stars coming from a region next to the bulge/bar intersection. Radial migration 
of the stars in the disk have been suggested by Haywood [25], who estimated upper val- 
ues for the migration rate from 1.5 to 3.7 kpc/Gyr, which agree with the values estimated 
in ll26l for the radial wandering due to the scattering mechanism assumed by Sellwood and 
Binney |fT9l. High resolution cosmological simulations I27J confirm that such scattering 
mechanism determines a significant migration in the stellar disk. However these simula- 
tions ignore the important effects of bars found by Minchev et al. ||20l[2Tl. Radial migration 
of stars (and gas) could have important implications for the interpretation of key observa- 
tional constraints, such as the age-metallicity relationship or the metallicity gradients, since 
old stars that formed at small galactocentric radii from enriched gas or young metal-poor 
stars at large radii are enabled to appear in a Solar-neighborhood sample. Due to the lack of 
detailed information on the process driving radial mixing, models of the Galactic chemical 
evolution have evaluated past history of the solar neighborhood and the formation and evo- 
lution of the abundance gradients assuming that radial mixing did not played an important 
role |[28l|29l[30l[3Tl[32l. Recently, Schoenrich and Binney |33 | explored the consequences 
of mass exchanges between annuli by taking into account the effect of resonant scattering of 
stars described before. This approach appears to be successful to replicate many properties 
of the thick disk in the Solar neighborhood without requiring any merger or tidal event |[34l . 
However, again in this case, the strong mixing mechanisms driven by bar resonances were 
not taken into account, casting thus doubts on some of their conclusions. 

In order to include the effect of radial migration in chemical evolution models and to 
gain a global (chemical + kinematical) vision of the processes at play in the galactic disks, 
many dynamical aspects need to be further investigated and in particular the role of the bar, 
that is the main non-axisymmetric component in disk galaxies. This is what we start to do 
in the present paper. In particular, we consider a marginally- stable disk which develops a 
central bar and spiral arms with the aim of quantitatively estimating the time and length 
scales of star diffusion in the radial direction. The final goal is to investigate how these 
characteristic scales evolve in time and how they depend on the activity of the bar. We 
present in detail the methods used for calculating the diffusion coefficient and time-scale. 
We will apply these methods to models which include the halo component and disks with 
different grade of stability in a forthcoming paper. 

The paper is organised as follows. In section 2 we describe the simulation and the 
relevant parameters. In section 3 we solve the diffusion equation in axisymmetric systems, 
we define the diffusion coefficient, the diffusion time-scale and the diffusion length-scale, 
and we describe how we estimate these quantities from the simulation results. In section 4 
we present our results and in section 5 we summarise our findings. 
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Figure 1 : Left: Toomre-parameter (solid blue line) and Araki (cTz/crr) parameter (dashed red line, labels 
on the right y-axis) at the final time. Middle: rotation curve. Right bar-strength time evolution. 



2 A/^-body simulation 

We have run self-consistent A/^-body simulations starting from a bar-unstable axisymmetric 
model. The initial configuration is marginally stable, i.e. the initial value of the Toomre 
parameter Qt = crrK/(336Gl.) is Qt ~ 1, where cr^ is the radial velocity dispersion of 
the disk component, G is the gravitational constant, E is the disk surface density, k is the 
epicycle frequency defined by - RdQ?/dR + where Q is the circular frequency 
related to the global gravitational potential 0(/^, z, t) in the disk plane z = by = 
(l/R)dO/dR. 

The initial mass distribution in our simulations corresponds to a superposition of a pair 
of axisymmetric Miyamoto-Nagai disks of mass Mb, Md, horizontal scales Ab + B,Ad + B, 
and identical scale-height B, 

^mn(R,z) = I (1) 

i=B,D ^r2 + + V52~r?)2 

The first component represents the bulge (B), while the second the disk (D) (35). The 
parameters have been set to Ab - 0.07 kpc, A^ = 1.5 kpc, B - 0.5 kpc, Mb/Md = 3/17. 
The initial particle positions and velocities are found by a pseudo-random draw following 
the density law corresponding to eq. ([T]), truncated to a spheroid of semi-axes R = 30 kpc, 
z = 10 kpc. The number of particles in the disk-bulge component is N = 4 ■ 10^. We 
then solve for the first velocity moment equations (Jeans equations) to find an approximate 
local equilibrium. The resulting distribution is then relaxed for a couple of rotations until 
ripples spreading through the disk from the center disappeared. We used this as the initial 
condition for the A/^-body simulations, performed with the Gadget-2 code |[36ll37l . 

The initial value of the radial and vertical velocity dispersions at two scale lengths 
from the center is cTr cr^ 20 km/s. The disk is marginally stable and spiral waves 
develop with the main global effect of heating in the radial direction. The Toomre and 
Araki parameters at the final time are shown in the left panel of fig. [T] solid blue and 
dashed red lines, respectively. The rotation curve at the final time is shown in the middle 
panel of fig. [1} 

We classify stellar orbits into three dynamical categories |[38l [35l: 1) the bar orbits 
and the disk orbits with the Jacobi integral H - E - QpL^ smaller than the value at the 
Lagrangian points Li 2, H < //(Li 2), where E is the total energy and is the z-component 
of the angular momentum. The particle separation in the bar or disk component can easily 
be done since the bar orbits have typically smaller values of and E. 2) The hot orbits 
with//>//(Li,2). 

The bar pattern speed Qp = dO/dt, where 6 is the azimuthal angle of the bar major 
axis (in the inertial frame) calculated by diagonalising the moment of inertia tensor of the 
bar particles, is 40 km s"^ kpc"^ We find that the pattern speed slowly decreases in time 
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Figure 2: Density maps at f ~ 550 Myr. Left: face-on view. Right: edge-on view. 



with a rate of a few km/s/kpc/Gyr in agreement with f39ll4]. The corresponding corotation 
radius Rc, obtained by the intersection of with the circular frequency, = ^{Rc, t), 
increases in time (it ranges from Rc -2 io A kpc in our simulation). 

In order to understand the role played by the central bar on the distribution of stars in 
the disk, we follow the evolution of the bar strength in time. The amplitude Cm of the mode 
m in the density distribution is defined as 



Cm — 



^ exp(/ mOj) 



(2) 



where j labels the stars. The bar's strength is defined as the mode Cm=2 when the stars 
are restricted to the bar component. We normalise C2 with respect to the number of stars 
in the bar component, Cq. This quantity is shown in the right panel of fig. [T] Since we 
do not include the gas component in our model, we are not able to follow its evolution for 
times longer than few Gyr. After this typical time-scale, the bar amplitude saturates and the 
system reaches a quasi-steady state. Since the inclusion of the gas component necessarily 
requires the introduction of other less controlled parameters, such as the cooling rate or star 
formation, we prefer to limit the integration time and to study the physical process at play 
within few Gyr. 

The face-on and edge-on views of the density distribution are shown in fig. [2] at time 
t ~ 600 Myr, that is just after the moment when the bar's strength reaches its maximum. 
Both bar and spiral arms develop, since the disk is sufficiently cold. In the external regions, 
the dominant pattern has m-\. 



3 Diffusion equation in axisymmetric systems 

We model the flow along the radial direction in the galactic plane by introducing the func- 
tion F{R, t) which satisfies the spatial diffusion equation 

dtF=^dR(RdRF), (3) 

K 

where D is the diffusion coefficient. The general solution of eq. ([3]) which is non-singular 
at = is given by 

F(R, t) =1 A{s) e'^"^ Jo(sR) s ds , (4) 
Jo 
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Figure 3: Two distributions with D = I and centered in Rq = 0.5 and 2 (solid lines) evolve in time 
(dashed and dotted lines, respectively), as described by eq. (|7]). 



where /o(-^) = Jo(-^) is the Bessel function of the first kind. The function A can be 
determined by taking the Hankel transform of F(R, 0) 



A(s) 



Jo 



F(R,0)Jo(sR)RdR. 



(5) 



Thus, by inserting eq. ^ into eq. ^ and assuming that the particles are initially localized 
at a certain radius Rq at time to = 0, F(R, 0) = FoRo6(R - Rq), we obtain 



F(R,t) = RlFo 

Jo 



s e'^"^ MsR) MsRo) ds . 



(6) 



Thus, the diffusion of this distribution can be expressed in term of elementary and Bessel 
functions iHOll . The time-evolution of an initial set of localized particles reads 



RiFo 
F(R, t) = exp 
2Dt ^ 



Rl + R^" 
4Dt 



(RRo\ 
\2Dtl ' 



(7) 



where Io(x) is the modified Bessel function of the first kind, which is finite at the origin, 
/q(0) = 1. In fig.[3]two initial distributions with D = I and centered in Rq = 0.5 and 2 (solid 
lines) evolve in time, as described by eq. ^ (dashed and dotted lines, respectively). At 
large radii, the distributions are essentially Gaussian, while at small radii they are strongly 
modified from the contributions of particles at the center of the cylinder. Eq. ^ describes 
the distribution of the radial positions of stars in the disk at the initial time ti which diffuse 
toward position Rq at time to, such that to - ti = At < To, where To is the diffusion time- 
scale or, equivalently, the distribution of stars which initially are in Ro at and then diffuse 
toward R with At < To- 

For R which goes to zero, eq. ^ reduces to 

RlFo 



2Dt 

For large values of the argument the modified Bessel function 1q{x) 
and thus F{R, f) reduces to 



(8) 



(2;rx)-i/2 exp(x) 



lim F{R, t) = 



R^oo 



^lAnDtR 



exp 



ADt 



yfR 



N(m,o-) 



(9) 




Figure 4: Left: radial dispersion of particles initially at (RqJo), for Rq = 2,4, . . . , 12 kpc and to = 
700, 1150, 1650 Myr, as a function of time. The growth rates in the vicinity of the initial time to give a 
first estimate of the diffusion coefficient from the relation (AR)^ - q-^ - 2Dt (dashed red lines). Right: 
distribution of particles at radius Ro = 2, 4, 6, 8 kpc and different initial times (dark blue), distribution 
of the same particles at time t, such that t - to < (light blue) and fit calculated by the nonlinear 
lest-square method (red line). 



where 7V(yu, cr) is the Gaussian distribution with mean value jd = (R) = Ro and standard 
deviation cr = ^/2Dt. 

As we will describe in the following subsection, the A/^-body simulation presented in 
section 2 allows us to estimate the diffusion coefficient D and the diffusion time-scale To = 
\to - ti\ in marginally stable disks. 



3.1 Calculation of diffusion time-scale and diffusion coefficient 

The diffusion coefficient D in eqs. Q-Q can be regarded as an instantaneous coefficient 
which depends on the position at which particles are initially localised and on the diffusion 
time. Modeling of stellar migration as a diffusion process is valid only for time intervals 
less than the diffusion time-scale. At < To, and for time intervals larger than the typical 
chaotic time-scale. Thus, the first issue is to estimate this diffusion time-scale from the 
simulation results. In the right panel of fig.|3] we show the radial dispersion (AR^) of stars 
initially at (Ro, to) as a function of time, for different values of (Ro, to). Ait = to the radial 
dispersion is null since the particles are all localised at R = Ro. As time evolves, (AR^) 
grows with a linear rate and then it saturates after a diffusion time-scale T^. The diffusion 
time scale turns out to be of the order of the rotation period, ~ T^ot = 27i/Q(R, t). 

From the relation (AR^) = cr^ = 2Dt, we calculate a first estimate of the diffusion 
coefficient D as the growth rates in the vicinity of the initial time to (red dashed lines in the 
left panel of fig. [3]) and we use these values as initial guess of the nonlinear least-square 
method, which minimizes the difference between the numerical results and the general 
solution of the diffusion equation described by eq. ^ for times less than T^. The maximum 
relative error on the estimated parameter D is AD/D ~ 0.06 near the central region. The 
fitting function that we obtain by this procedure reproduces reasonably well the overall 
distribution of particles at times t > to (or t < to), with At = \t - to\ < T^, as can be seen in 
the right panel of fig. |4j 
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Figure 5: Left: contour map of the diffusion coefficient D. Right: L^-values of the particles at time 
2.2 Gyr within the radial ranges R = (2.5 ± 0.5) kpc (bin 1, blue line), R = (3.0± 0.5) kpc (bin 2, red 
line) and R = (8.0 + 0.5) kpc (bin 3, black line). 



4 Results 

The methods described in the previous section allowed us to calculate the diffusion coef- 
ficient D(R, t), which is shown in the contour map in the left panel of of fig. [5] The bar's 
corotation radius is also shown as a dashed white line. We can see that the diffusion coeffi- 
cient is not constant in time nor in radius. If the disk is not too hot, D has the largest values 
D ~ 0.1 kpc^Myr"^ near the corotation radius of the bar (which increases in time in our 
simulation), where the density is strongly perturbed by a m = 2 pattern created by the bar 
and the transient spiral arms, and in the external regions > 8 kpc, where the density is 
modulated by a m = 1 pattern. The stars respond collectively to these modulations and the 
process of migration corresponds to a diffusion in an axisymmetric system. 

From D we can estimate the radial dispersion cr in a rotation period (which is of the 
order of the diffusion time-scale), cr = ^^IDTrot- If the disk is sufficiently cold, the radial 
dispersion is high near the corotation region of the bar, where it can recurrently assume 
values of the order of cr ~ 6 kpc. At an intermediate radius, such as = 6 kpc, the radial 
dispersion is of the order of cr ~ 3 kpc and it increases in the external less dense regions 
where the pattern m = \ dominates. In this external region, the radial migration is very 
high and the number of stars which stay in a local volume of radius d ^ 100 pc around, for 
example, = 8 kpc is very low for times less than the diffusion time, since d ^ cr ^ 5 kpc. 

Three different bins of particles located respectively 'mR = (2.5 + 0.5) kpc, = (3.0 + 
0.5) kpc, R = (8.0 ± 0.5) kpc ait = 2.2 Gyr are followed in time and their evolution history 
is shown in fig.|6] In the first bin, particles are inside the bar and they remain there all along 
the evolution. The presence of a m = 1 mode in the central bar region can be observed 
as periodical modulations in the density. These particles have values of ~ (see the 
right panel of fig. [5} blue line labeled as 'bin 1') and large negative energies (see the left 
panel of fig.|7]). Particles in the second bin centered inR = (3.0 + 0.5) kpc at / = 2.2 Gyr 
belong to two different types of orbits: one family can migrate only inside, the other can go 
outside the bar, in the disk. Large values of the diffusion coefficient D near the corotation 
region (see the left panel of fig. [5]) are related to this superposition of two families of bar 
orbits. The corresponding values of and E are shown, respectively, in the red line in fig. [5] 
(left panel, line labeled as 'bin 2') and in the middle panel of fig.[7| where it can be easily 
seen that the two families have two different ranges of energies: large negative energies 
correspond to bar particles and small negative energies correspond to disk particles. In 
the third bin, particles belong to the disk component: they can span all the region outside 
corotation, Rc < R < 12 kpc, and they come mainly from the corotation region at the time 
when the bar strength had a maximum (i.e. at / ~ 350 Myr, cf. fig.[T]l). These disk particles 
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Figure 6: Radial migration of particles which at time t = 2.2 Gyr are inR = (2.5 + 0.5) kpc (left panel), 
R = (3.0 + 0.5) kpc (middle panel) and R = (8.0 + 0.5) kpc (right panel). 
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Figure 7: Energies of particles which at time t = 2.2 Gyr are in R = (2.5 + 0.5) kpc (left panel), 
R = (3.0 + 0.5) kpc (middle panel) and R = (8.0 + 0.5) kpc (right panel). 

have small negative energies (see the right panel of fig. |7]) and large values of (see the 
black line labeled as 'bin 3' in the right panel of fig. [5]). 

5 Conclusions 

Modeling the migration of stars in marginally stable disks as a diffusion process in the 
radial direction is a powerful tool which allows us to estimate quantitatively two crucial 
parameters, the diffusion coefficient and the diffusion time-scale. It is important to note 
that such diffusion model is only valid for describing the system for times smaller than the 
diffusion time-scale, which is of the order of the rotation period, and for times larger than 
the chaotic time-scale, after which orbits are no more time-reversible. 

The calculations of both the diffusion coefficient and the diffusion time-scale give us a 
quantitative measure of the migration process in the disk. Another advantage of studying 
the diffusion of stars in real space, rather than in velocity space, is that it can be easily 
related to the evolution of chemical elements, thus representing an interesting tool to be 
implemented in chemical evolution codes. 

We have found that the diffusion time-scale is of the order of one rotation period and 
that the diffusion coefficient D depends on the evolution history of the disk and on the radial 
position. Larger values D are found near the corotation region, which evolves in time, and 
in the external region, where asynmietric patterns develop. Marginally stable disks, with 
2r ~ 1, have two different families of bar orbits with different values of angular momentum 
Lz and energy E, which determine a large diffusion in the corotation region. 

The application of the method described here to models which include the halo com- 
ponent and have disks with different initial Toomre-parameter Qt will be presented in a 
forthcoming paper. 
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